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Abstract. The main reason for the generation of tsunamis is the deformation of the 
bottom of the ocean caused by an underwater earthquake. Usually, only the vertical 
bottom motion is taken into account while the horizontal co-seismic displacements are 
neglected in the absence of landslides. In the present study we propose a methodology 
based on the well-known Okada solution to reconstruct in more details all components of 
the bottom coseismic displacements. Then, the sea-bed motion is coupled with a three- 
dimensional weakly nonlinear water wave solver which allows us to simulate a tsunami 
wave generation. We pay special attention to the evolution of kinetic and potential energies 
of the resulting wave while the contribution of the horizontal displacements into wave 
energy balance is also quantified. Such contribution of horizontal displacements to the 
tsunami generation has not been discussed before, and it is different from the existing 
approaches. The methods proposed in this study are illustrated on the July 17, 2006 Java 
tsunami and some more recent events. 
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1. Introduction 

During the years 2004 to 2006 several interesting tsunami events took place in the Indian 
Ocean. In December 26, 2004, the Great Sumatra-Andaman earthquake (M w = 9.1, cf. 
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C. Amnion et al. (2005), [1]) generated the devastating Indian Ocean tsunami refered to 
in the literature as the Tsunami Boxing Day 2004 (cf. C. Synolakis & E. Bernard (2006), 
[83]). The local tsunami run-ups from this event exceeded the height of 34 m at Lhoknga 
in the western Aceh Province. This and other observations led many researchers to ask 
whether this tsunami was unusualy large for this specific earthquake size. Later it was 
shown by E.L. Geist et al. (2006) [32] that the Great Sumatra-Andaman earthquake is 
very similar in terms of local tsunami magnitude to past events of the same size. For 
example, the 1964 Great Alaska earthquake (M w = 9.2, cf. H. Kanamori (1970) [49]) 
demonstrates a similar scaling. 

On March 28, 2005 another earthquake occured approximately 110 km to the SE from 
the Great Sumatra-Andaman earthquake's epicenter. The magnitude of this earthquake 
was estimated to be M w = 8.6 ~ 8.7 (cf. [57, 8, 93]). This event triggered a massive 
evacuation in the surrounding Indian Ocean countries. However, the March 28, 2005 
Northern Sumatra earthquake failed to generate a significant tsunami event. The survey 
teams reported maximum tsunami run-up of 4 m. [32]. This event can be compared with 
the 1957 Aleutian earthquake (M w = 8.6, cf. J.M. Johnson and K. Satake (1993) [46]) 
which produced a maximum tsunami run-up of 15 m (see J.F. Lander (1996) [55]). The 
deficiency of the March 2005 tsunami is related to the slip concentration in the down-dip 
part of the rupture zone and to the fact that a substantial part of the vertical displacement 
occurred in shallow waters or on the substance of the ground [32]. 

On the other hand, the smaller July 16, 2006 Java earthquake (M w = 7.8, cf. C.J. Am- 
nion et al. (2006), [2]) generated an unexpectedly destructive tsunami wave which affected 
over 300 km of south Java coastline and killed more than 600 people, [26]. Field mea- 
surements of run-up distributions range uniformly from 5 to 7 m in most inundated areas. 
However, unexpectedly high run-up values were reported at Nusakambangan Island exceed- 
ing the value of 20 m. This tsunami focusing effect could be seemingly ascribed to local 
site effects and/or to a local submarine landslide/slump. July 16, 2006 Java tsunami can 
be compared to a similar event occured on June 2, 1994 at the East Java (M w = 7.6) (see 
Y. Tsuji et al. (1995), [88]). This 1994 Java tsunami produced more than 200 casualties 
with local run-up at Rajekwesi slightly exceeding 13 m. 

All these examples of recent and historical tsunami events show that there is a big variety 
of the local tsunami heights and run-up values with respect to the earthquake magnitude 
M w . It is obvious that the seismic moment M of underwater shallow earthquakes is 
adequate to estimate far- field tsunami amplitudes (see also [68]). However, tsunami wave 
energy reflects better the local tsunami severity, while the specific run-up distribution 
depends on bathymetric propagation paths and other site-specific effects. One of the 
central challenges in the tsunami science is to rapidly assess a local tsunami severity from 
very first rough earthquake estimations. In the current state of our knowledge false alarms 
for local tsunamis appear to be unavoidable. The tsunami generation modeling attempts 
to improve our understanding of tsunami behaviour in the vicinity of the genesis region. 

Tsunami generation modeling was initiated in the early sixties by the prominent work 
of K. Kajiura, [47], who proposed the use of the static vertical sea bed displacement 
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for the initial condition of the free surface elevation. Classically, the celebrated Okada, 
[66, 67], and sometimes Mansinha & Smylie 1 [59, 60] or Gusiakov,[33], solutions are used 
to compute the coseismic sea bed displacements. This approach is still widely used by the 
tsunami wave modeling community. However, some progress has been recently made in 
this direction, [65, 19, 18, 21, 74, 76, 25]. 

There is a consensus on the importance of the horizontal motion in landslide generated 
tsunamis, [37, 95, 4, 69]. However, for tsunami waves caused by underwater earthquakes, 
the horizontal displacements are often disregarded by the tsunami wave modeling com- 
munity. We can quote here a few publications devoted to tsunami waves such as one by 
D. Stevenson (2005), [81]: 

"Although horizontal displacements are often larger, they are unimportant 
for tsunami generation except to the extent that the sloping ocean floor also 
forces a vertical displacement of the water column." 

Or another one by G.A. Ichinose et al. (2000) [41]: 

"The initial lake level values were specified by assuming that the lake surface 
instantaneously conformed to the vertical displacement of the lake bottom 
while horizontal velocities were set to zero. The effect of horizontal defor- 
mation on the initial condition is neglected here and left for future work." 

The authors of this article also tended to neglect the horizontal displacements field in 
previous tsunami generation studies, [18, 64]. Perhaps, this situation can be ascribed to 
the work of E. Berg (1970), [7], who showed in the case of the 1964 Alaska's earthquake 
(M w = 9.2) that the input into the potential energy from the horizontal motion is less than 
1.5% of that from the vertical movement. 

The attitude to horizontal displacements changed after the prominent work by Y. Tan- 
ioka and K. Satake (1996), [84]. According to their computations for the 1994 Java earth- 
quake, the inclusion of horizontal displacements contribution results in 43% increase in 
maximum initial vertical ground displacement and an increase in the wave amplitudes at 
the shoreline by 30%. These striking results incited researchers to reconsider the role of 
the horizontal sea bed motion. Hereafter, various tsunami generation techniques which 
involve only the vertical displacement field are referred to as the incomplete scenario. On 
the contrary, the complete tsunami generation takes also into account the horizontal dis- 
placements field. The term complete in our study refers to a purely kinematic interaction 
between the moving bottom and the ocean layer in contrast to the horizontal impulse- 
momentum transfer mechanism [91] applied for the first time to tsunami wave generation 
in [78]. 

The approach proposed by Tanioka & Satake (1996) is based on a simple physical con- 
sideration: horizontal displacements of a sloping bottom will produce some amount of the 
vertical motion depending on the slope magnitude. Assuming that the slope is small, this 
idea can be expressed mathematically by applying once the first order Taylor expansion: 

u h := Uih x + u 2 h v , 



In fact, Mansinha & Smylie solution is a particular case of the more general Okada solution. 
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where uxp are the horizontal components of the displacements field and h(x, y) is the 
bathymetry function. The indices x and y denote the spatial partial derivatives. Tanioka 
& Satake (1996) referred to the quantity Uh as the vertical displacement of water due to the 
horizontal movement of the slope. Then, the quantity Uh is simply added as a correction to 
the vertical displacements component us(x,y). To compute the co-seismic displacements 
u i,2,3 they used the celebrated Okada solution [80, 66, 67] applied to a single rectangular 
fault segment. In the present study we apply the same solution to multiple sub-faults 
to achieve higher resolution (see also [24]). More recently, an original approach based on 
the impulse-momentum principle, [91], has been proposed in [78]. Their genesis theory 
is based on the idea of the momentum exchange between the slipped continental slope 
and the water column due to the earthquake. The horizontal bottom motion contributes 
directly to the ocean kinetic energy as any moving object in water transfers momentum 
to the fluid. This idea is similar to the classical added mass concept. According to the 
authors of [78], the continental slope is not an exception from this principle. This implies 
some special treatment of the water parcels in the vicinity of moving boundaries. For more 
details we refer to the original publication [78]. The horizontal impulse approach leads to 
predictions which are somehow different from the present and several previous studies. The 
methodology described in this article is closer to the ideas of Tanioka & Satake (1996) [84]. 
Namely, we reconstruct in more details the kinematics of bottom displacements in vertical 
and horizontal directions, including their evolution in time. We underline that most of 
tsunami generation studies, including [84], assume the bottom deformation process to be 
instantaneous, that in most of the cases provides a fairly good approximation. Then, the 
bottom motion generates some perturbations on the free surface of the water layer. This 
perturbation form long waves, propagating across the oceans under the gravity force and 
commonly referred to as tsunami waves due to their destructive potential fully realized 
only in coastal regions, [47, 89, 16, 83]. 

In the current study we shed some light into the energy transfer process from an under- 
water earthquake to the implied tsunami wave (in the spirit of the study by D. Dutykh 
& F. Dias (2009), [20]), while taking into account and quantifying the horizontal displace- 
ments contribution into tsunami energy balance. We focus only on the generation stage 
since the propagation phase and run-up techniques are well understood nowadays at least 
in the context of Nonlinear Shallow Water equations, [42, 14, 53, 25, 23]. 

The present study is organized as follows. In Section 2 we present mathematical models 
used in this study. Specifically, in Sections 2.1 and 2.2 a description of the bottom motion 
and the fluid layer solution is presented respectively. Some rationale on tsunami wave 
energy computations is discussed in Section 2.2.1. Numerical results are presented in 
Section 3. Finally, some important conclusions of this study are outlined in Section 4. In 
the Appendix of this paper we present the results concerning the tsumami generation of 
two recent events. 
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2. Mathematical models 

Tsunami waves are caused by a huge and rapid motion of the seafloor due to an under- 
water earthquake over broad areas in comparison to water depth. There are some other 
mechanisms of tsunami genesis such as underwater landslides, for example. However, in 
this study we focus on the purely seismic mechanism which occurs most frequently in 
nature. 

Hydrodynamics and seismology are only weakly coupled in the tsunami generation prob- 
lem. Namely, the released seismic energy is partially transmitted to the ocean through the 
sea bed deformation while the ocean has obviously no influence onto the rupturing process. 
Consequently, our problem is reduced to two relatively distinct questions: 

(1) Reconstruct the sea bed deformation h = h(x,t) caused by the seismic event under 
consideration 

(2) Compute the resulting free surface motion 

The answers to these questions are analyzed in Sections 2.1 and 2.2 respectively. 

2.1. Dynamic bottom displacements reconstruction. Traditionally, the free surface 
initial condition for various tsunami propagation codes (see [85, 42, 43]) is assumed to be 
identical to the static vertical deformation of the ocean bottom, [48]. This assumption is 
classically justified by the three following reasons: 

(1) Tsunamis are long waves. In this regime the vertical acceleration is neglected with 
respect to the gravity force 

(2) The sea bed deformation is assumed to be instantaneous. It is based on the com- 
parison of gravity wave speed (200 m/s for water depth of 4 km) and the seismic 
wave celerity (~ 3600 m/s) 

(3) The effect of horizontal bottom motion is negligible for tsunami generation since 
the bathymetry has in general mild slope (~ 10%), [7] 

It is worth to note that nonhydrostatic effects as well as finite time source duration have 
been modeled in several recent studies [87, 19, 28, 52, 21, 29]. Some attempts have also 
been made to take into account the horizontal displacements contribution, [84, 32, 78]. 

In the present study we relax all three classical assumptions. The bottom deformation 
is reconstructed using the finite fault solution as it was suggested in our previous study 
[24]. However, we extend the previous construction to take into account the horizontal 
displacements contribution as well. The finite fault solution is based on the multi-fault 
representation of the rupture, [5, 45]. The rupture complexity is reconstructed using the 
inversion of seismic data. Fault's surface is parametrized by multiple segments with variable 
local slip, rake angle, rise time and rupture velocity. The inversion is performed in an 
appropriate wavelet transform space. The objective function is a weighted sum of L%, Li 
norms and some correlative functions. With this approach seismologists are able to recover 
rupture slip details, [5, 45]. This available seismic information is exploited hereafter to 
compute the sea bed displacements produced by an underwater earthquake with better 
geophysical resolution. Several other multiple segments sources have been used in [31, 94, 
61]. 
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Figure 1. Surface projection of the fault's plane and the ETOPOl bathy- 
metric map of the region under consideration. The symbol * indicates the 
hypocenter's location at (107.345°, —9.295°). The local Cartesian coordinate 
system is centered at the point (108°, —10°). This region is located between 
(106°, —8°) and (110°, —12°). The colorbar indicates the water depth in me- 
ters below the still water level (z = 0). In the region under consideration 
the depth varies from 20 to 7100 meters. 

Remark 1. In a few studies an attempt has been made to reconstruct the seismic source 
from tsunami tide gauge records, [71, 72, 97, 30, 40, 75]. This approach seems to be very 
promising and in future a joint combination of seismic and hydrodynamic inversions should 
be used for the successful reconstruction of appropriate initial conditions. 

Let us describe the way of how the sea bed displacements are reconstructed. In this 
reconstruction procedure we follow the great lines of our previous study [24], while adding 
new ingredients concerning the horizontal displacements field reconstruction. We illustrate 
the proposed approach in the case of the July 17, 2006 Java tsunami for which the finite 
fault solution is available, cf. [44, 70]. It was a relatively slow earthquake and thus, 
atypical. However, we assume that the slow rupturing process was well resolved by the 
finite fault inversion algorithm. 

The fault is considered to be the rectangle with vertices located at (109.20508° (Lon), 
-10.37387° (Lat), 6.24795 km (Depth)), (106.50434°, -9.45925°, 6.24795 km), (106.72382°, 
-8.82807°, 19.79951 km), (109.42455°, -9.74269°, 19.79951 km) (see Figure 1). The fault's 
plane is conventionally divided into N x = 21 subfaults along strike and N y = 7 subfaults 
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P-wave celerity c p , m/s 


6000 


S-wave celerity c s , m/s 


3400 


Crust density p, kg/m 3 


2700 


Dip angle, 5 


10.35° 


Slip angle (CW from N) 


288.94° 



Table 1. Geophysical parameters used to model elastic properties of the 
subduction zone in the region of Java. 

down the dip angle, leading to the total number of N x x N y = 147 equal segments. Pa- 
rameters such as the subfault location (x c , y c ), the depth d iy the slip Ui and the rake angle 
4>i for each segment can be found in [44] (see also Appendix II, [24]). The elastic constants 
and parameters such as dip and slip angles, which are common to all subfaults, are given 
in Table 1. We underline that the slip angle is measured conventionally in the counter- 
clockwise direction from the North. The relations between the elastic wave celerities c p , c s 
and Lame coefficients A, /i used in Okada's solution are classical and can also be found in 
Appendix III, [24]. 

One of the main ingredients in our construction is the so-called Okada solution, [66, 67], 
which is used in the case of an active fault of small or intermediate size. The success of this 
solution may be ascribed to the closed-form analytical expressions which can be effectively 
used for various modeling purposes involving co-seismic deformation. 

Remark 2. The celebrated Okada solution, [66, 67], is based on two main ingredients 
- the dislocation theory of Volterra [92, 80] and Mindlin's fundamental solution for an 
elastic half-space, [63]. Particular cases of this solution were known before Okada's work, 
for example the well-known Mansinha and Smylie's solution, [59, 60]. Usually all these 
particular cases differ by the choice of the dislocation and the Burger's vector orientation, 
[73]. We recall the basic assumptions behind this solution: 

• Fault is immersed into the linear homogeneous and isotropic half-space 

• Fault is a Volterra's type dislocation 

• Dislocation has a rectangular shape 

For more information on Okada's solution we refer to [19, 16, 18] and the references therein. 

The trace of the Okada's solution at the sea bottom (substituting z = in the geophysical 
coordinate system) for each subfault will be denoted o!f\x; 5, A, /i, . . .), where 5 is the dip 
angle, A, /i are the Lame coefficients and dots denote the dependence of the function Of' [x] 
on other 8 parameters, cf. [19]. The index i takes values from 1 to N x x N y and denotes the 
corresponding subfault segment, while the superscript j is equal to 1 or 2 for the horizontal 
displacements and to 3 for the vertical component. Hereafter we will adopt the short-hand 
notation 0^ (x) for the j th displacement component of the Okada's solution for the i th 
segment having in mind its dependence on various parameters. 

Taking into account the dynamic characteristics of the rupturing process, we make some 
further assumptions on the time dependence of the displacement fields. The finite fault 
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solution provides us with two additional parameters concerning the rupture kinematics - 
the rupture velocity v T and the rise time t r which are equal to 1.1 km/s and 8 s for July 17, 
2006 Java event respectively The epicenter is located at the point x e = (107.345°, —9.295°), 
cf. [44]. Given the origin x e , the rupture velocity v r and 2 th subfault location Xi, we define 
the subfault activation times ti needed for the rupture to achieve the corresponding segment 
i by the formulas: 

U = — — — , i — l,...,N x x N y . 

Consequently, for the sake of simplicity in our study we assume the rupture propagation 
along the fault to be isotropic and homogeneous. However, some more sophisticated ap- 
proaches including possible heterogeneities exist (see e.g. [13]). We will also follow the 
pioneering idea of J. Hammack, [36, 35], developed later in [87, 86, 19, 22, 52], where the 
maximum bottom deformation is achieved during some finite time (the so-called rise time) 
according to an appropriately chosen dynamic scenario. Various scenarios used in practice 
(instantaneous, linear, trigonometric, exponential, etc) can be found in [36, 50, 35, 22, 19]. 
In this study we adopt the trigonometric scenario which is given by the following formula: 

T(t) =H(t-t r ) + ~H(t)H(t r -t)(l-cos('Kt/t r )), (2.1) 

where %{t) is the Heaviside step function. This scenario has the advantage to have also the 
first derivative continuous at the activation time t = 0. However, for comparative purposes 
sometimes we will use also the so-called exponential scenario [50]: 

T e {t)=H{t)(l-e- at ), «:=^P- 

For illustrative purposes both dynamic scenarios are represented in Figure 2. 

We sum up together all the ingredients proposed above to reconstruct dynamic displace- 
ments field u = (ui,U2,us) at the sea bottom: 



Uj{X, 



t)= £ T(t-U)0?\x). 



i=l 



Remark 3. We would like to underline here the asymptotic behaviour of the sea bed dis- 
placements. By definition of the trigonometric scenario (2.1) we have lim T(t) = 1. 

Consequently, the sea bed deformation will attain fast its state which consists of the linear 
superposition of subfaults contributions: 



i=l 



Finally, we can predict the sea bed motion by taking into account horizontal and vertical 
displacements: 

h(x, t) = h (x - Mi >2 (x , t)) - u 3 (x, t), (2.2) 
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Dynamic scenarios 




t, s 

Figure 2. Trigonometric and exponential dynamic scenarios for t r = 1 s 
(see J. Hammack (1973), [35]). 

where ho(x) is a function which interpolates 2 the static bathymetry profile given e.g. by 
the ETOPOl database (see Figures 1 and 4). 

Remark 4. In some studies where horizontal displacements were taken into account (cf. 
[84, 3, 78, 6}), the first order Taylor expansion was permanently applied to the bathymetry 
representation formula (2.2) to give: 

h(x, t) w h (x) - ui, 2 (x, t) ■ Vgho(x) - u 3 (x, t) (2.3) 

We prefer not to follow this approximation and to use the exact formula (2.2) since it 
is valid for all slopes (see Figure 4 f or Java bathymetry). Another difficulty lies in the 
estimation of the bathymetry gradient Vgho(x) required by Taylor's formula (2.3). The 
application of any finite differences scheme to the measured h (x) leads to an ill-posed 
problem. Consequently, one needs to apply extensive smoothing procedures to the raw data 
ho(x) which induces an additional loss in accuracy. 

In the present study we do not completely avoid the computation of the static bathymetry 
gradient Vgho(x) since the kinematic bottom boundary condition (2.7) involves the time 
derivative of the bathymetry function: 

d t h = -V$h (x) ■ d t u 1)2 (x,t) - d t u 3 (x,t). 



In our numerical simulations presented below we use the MATLAB TriScatteredlnterp class to 
interpolate the static bathymetry values given by ETOPOl database. 
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However, the last formula is exact and it is obtained by a straightforward application of the 
chain differentiation rule. 

Another possibility could be to consider static horizontal displacements u\^{x) thus keep- 
ing dynamics only in the vertical component u^{x,t). However, we do not choose this 
option in this work. 

In the next section we will present our approach in coupling this dynamic deformation 
with the hydrodynamic problem to predict waves induced on the free surface of the ocean. 



2.2. The water wave problem with moving bottom. We consider the incompressible 
flow of an ideal fluid with constant density p in the domain O C M 2 . The horizontal 
independent variables will be denoted by x = (x, y) and the vertical one by z. The origin 
of the Cartesian coordinate system is traditionally chosen such that the surface z = 
corresponds to the still water level. The fluid domain is bounded below by the bottom 
z = —h(x,t) and above by the free surface z = rj(x,t). Usually we assume that the total 
depth H(x,t) := h(x,t) + r](x,t) remains positive H(x,t) > h > under the system 
dynamics V£ G [0, T]. The sketch of the physical domain is shown in Figure 3. 

Remark 5. Classically in water wave modeling, we make the assumption that the free 
surface is a graph z = rj(x,t) of a single-valued function. It means in practice that we 
exclude some interesting phenomena, (e.g. wave breaking phenomena) which are out of the 
scope of this modeling paradigm. 

The governing equations of the classical water wave problem are the following, [54, 82, 
62, 96]: 



--V 2 <j) + d 2 zz <J) = 0, (x,z)enx[-h,ri\, (2.4) 

d t r] + V0 • V?7 - d z (j) = 0, z = rj(x, t), (2.5) 

<9 t + ±|V0| 2 + f(<9 z 0) 2 + ^ = 0, z = V (x,t) } (2.6) 

d t h + V0 • Vh + d z (f> = 0, z = -h(x, t), (2.7) 

where is the velocity potential, g the acceleration due to gravity force and V = (d x , d y ) 
denotes the gradient operator in horizontal Cartesian coordinates. The fluid incompressibil- 
ity and flow irrotationality assumptions lead to the Laplace equation (2.4) for the velocity 
potential 4>(x, z, t). 

The main difficulty of the water wave problem lies on the boundary conditions. Equations 
(2.5) and (2.7) express the free surface and bottom impermeability, while the Bernoulli 
condition (2.6) expresses the free surface isobarity respectively. 

Function h(x, t) represents the ocean's bathymetry (depth below the still water level, 
see Figure 3) and is assumed to be known. The dependence on time is included in order to 
take into account the bottom motion during an underwater earthquake [16, 19, 22, 52, 18]. 

For the exposition below we will need also to compute unitary exterior normals to the 
fluid domain. The normals at the free surface Sf and at the bottom Sf, are given by the 
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n (x.y.t) 



h(x,y,t) 



Figure 3. Sketch of the physical domain. 



following expressions respectively: 

1 



n 



f 



-Vrj 
1 



n b 



-Vh 
-1 



In 1968 V. Zakharov proposed a different formulation of the water wave problem based 
on the trace of the velocity potential at the free surface [98]: 

<p(x,t) := (j)(x,r](x,t),t). 

This variable plays a role of the generalized momentum in the Hamiltonian description of 
water waves [98, 15]. The second canonical variable is the free surface elevation r\. 

Another important ingredient is the normal velocity at the free surface v n which is 
defined as: 



fa,t) := VT+Wrt 



dhf 



Z=T) 



(2.8) 



z=r) 



Kinematic and dynamic boundary conditions (2.5), (2.6) at the free surface can be rewritten 
in terms of ip, v n and rj [12, 11, 27]: 



dtV ~ V v ((p) = 0, 

dt<p + 1| v v \ 2 + 97]- 2(1+l x v?? | 2) ]pM + • v v ] 2 = o. 



(2.9) 
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Here we introduced the so-called Dirichlet-to-Neumann operator (D2N) [10, 11] which maps 
the velocity potential at the free surface ip to the normal velocity v n : 



Z = T\ 



V 2 <P + d 2 zz <P = 0, (x,z)eQx[-h,r)}, 
</> =tp, z = rj, 



l + \Vh\ 2 -f- =d t h, z = -h. 

The name of this operator comes from the fact that it makes a correspondance between 
Dirichlet data 09 and Neumann data \/l + I V?]! 2 — — at the free surface. 

J Z=1J 

So, the water wave problem can be reduced to a system of two PDEs (2.9) governing 
the evolution of the canonical variables rj and (p. For the tsunami generation problem we 
approximate and we compute efficiently the D2N map T> v (ip) using the Weakly Nonlinear 
(WN) model described in [24]. This relies on the approximate solution of the 3D Laplace 
equation in a perturbed strip-like domain using the Fourier transform (0 := J-[ip], r\ = 

Vr,{(p) = 0\k\tanh(\k\H) + f sedx(\k\H) - ^^[ikip] ■ J 7 " 1 [ikfj] 

where k is the wavenumber and / is the bathymetry forcing term 

f(x,t) :=-d t h-V<P\ z= _ h -Vh. 

Several details on the time integration procedure can be also found in [24]. The resulting 
method is only weakly nonlinear and analogous at some point to the first order approxima- 
tion model proposed in [34]. As the hydrodynamical models, Tanioka and Satake (1996) 
[84] used the Linear Shallow Water equations [77], while Song et al. (2008) employed a 
hydrostatic ocean circulation model [79]. 

Recently, it was shown that a tsunami generation process is essentially linear, [52, 76]. 
However, our WN approach will take into account some nonlinear effects when they become 
important, for example when rapid changes in the bathymetry are present. This is possible 
when the generation region involves a wide range of depths from deep to shallow regions 
(see Figures 1 and 4). 

2.2.1. Tsunami wave energy. In this study we are more particularly interested in the evolu- 
tion of the generated wave energy [20]. In the case of free surface incompressible flows, the 
kinetic and potential energies, denoted by K and n respectively, are completely determined 
by the velocity field and the free surface elevation: 

v 

K{t) :=§/ J J |V0| 2 dxdz, H(t) := ^ j j rf dx. 

-h n n 

The definition of the kinetic energy K{t) involves an integral over the three dimensional 
physical domain Qx [—h,rj\. We can reduce the integral dimension using the fact that the 
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velocity potential is a harmonic function: 

|V0| 2 = V • (0V0) - (j) A(f) = V ■ (0V0). 

=0 

Consequently, the kinetic energy can be rewritten as follows: 

K{t)= 2 // ^^■ fldo = \ fj <®+% // Hhdx, 

s f +s b n 



(i) (ii) 

where denotes the trace of the velocity potential at the bottom 4>\ z= _ h (see D. Clamond 
& D. Dutykh (2012), [9]). In order to obtain the last equality we used the free surface and 
the bottom kinematic boundary conditions (2.6), (2.7). The first integral (I) is classical 
and represents the change of kinetic energy under the free surface motion while the second 
one (II) is the forcing term due to the bottom deformation. The total energy 3 is defined 
as the sum of kinetic and potential ones: 

E(t) := K(t) + II(t) = ^JJ tpV^ip) dx+^JJ^dthdx + ^-JJr, 2 dx. 

n si u 

Below we will compute the evolution of the kinetic, potential and total energies beneath 
moving bottom. 

3. Numerical results 

The proposed approach will be directly illustrated on the Java 2006 event. The July 
17, 2006 Java earthquake involved thrust faulting in the Java's trench and generated a 
tsunami wave that inundated the southern coast of Java, [2, 26]. The estimates of the size 
of the earthquake, [2], indicate a seismic moment of 6.7 x 10 20 N • m, which corresponds to 
the magnitude M w = 7.8. Later this estimation was refined to M w = 7.7, [44]. Like other 
events in this region, Java's event had an unusually low rupture speed of 1.0 - 1.5 km/s 
(we take the value of 1.1 km/s according to the finite fault solution [44]), and occurred 
near the up-dip edge of the subduction zone thrust fault. According to C. Ammon et al, 
[2] , most aftershocks involved normal faulting. The rupture propagated approximately 200 
km along the trench with an overall duration of approximately 185 s. The fault's surface 
projection along with ocean ETOPOl bathymetric map are shown in Figures 1 and 4. We 
note that Indian Ocean's depth of the region considered in this study varies between 7186 
and 20 meters in the shallowest regions which may imply local importance of nonlinear 
effects. 

Remark 6. We have to mention that the finite fault inversion for this earthquake was 
also performed by the Caltech team, [70]. The main differences with the USGS inversion 
consist on the employed dataset. To our knowledge, A. Ozgun Konca and his collaborators 
include also displacements measured with GPS-based techniques. Consequently, they came 



3 We note that the total energy is not conserved during the tsunami generation phase due to the forcing 
term (II) coming from the bottom kinematic boundary condition. 
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Figure 4. Side view of the bathymetry, cf. also Figure 1. 



Ocean water density, p, kg/m 3 


1027.0 


Gravity acceleration, g, m/s 2 


9.81 


Epicenter location (Lon, Lat) 


(107.345°, -9.295°) 


Rupture velocity, v r , km/s 


1.1 


Rise time, t , s 


8.0 


Number of Fourier modes in x 


256 


Number of Fourier modes in y 


256 



Table 2. Values of physical and numerical parameters used in simulations. 

to the conclusion that the July 17, Southern Java earthquake magnitude was M w = 7.9. 
The energy of a tsunami wave generated according to this solution will be discussed below. 

The numerical solutions presented below are given by the Weakly Nonlinear (WN) model. 
A uniform grid of 256 x 256 points 4 is used in all computations below. The time step At 
is chosen adaptively according to the RK(4,5) method proposed in [17]. The problem 
is integrated numerically during T = 255 s which is a sufficient time interval for the 
bottom to take its final shape (< 220 s) and of the resulting tsunami wave to enter into 



Since we use a pseudo-spectral method, the convergence is expected to be exponential and this number 
of harmonics should be sufficient to capture all scales important for phenomena that we consider here. 
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Figure 5. Contribution of the horizontal displacements field into the sea 
bed deformation expressed in percentage of the maximum vertical displace- 
ment. The grey color corresponds to negligible values while white and black 
zones show the most important contributions. The colorbar shows the per- 
centage of the signed difference going from -21% to +7%. 

the propagation stage. The values of various physical and numerical parameters used in 
simulations are given in Table 2. 

We begin our numerical investigations by quantifying the contribution of horizontal dis- 
placements into the sea bed deformation process. For this purpose we consider the differ- 
ence dh{x) between the deformed bottom under the action of only horizontal displacements 
in their steady state (t — > +oo) and the initial configuration: 

, ._ hp(x- Ui, 2 (x,t)\ t ^ +oQ ) - h (x) 
max \u 3 {x,t)\ 

x,t—t+oo 

In Figure 5 we present the quantitative effect of horizontal displacements relative to the 
maximum vertical displacement max 1 7/3 (if, t = +oo)|. The computations we performed 

X 

show that the maximum amplitude of the bottom variation due to the action of horizontal 
displacements reaches 21% of the maximum amplitude of the vertical displacement. In 
practice it means that locally (depending on the bathymetry shape and the slip distribu- 
tion) we cannot completely neglect the effect of horizontal motion. 

The next step consists in quantifying the impact of horizontal displacements onto free 
surface motion. We put six numerical wave gauges at the following locations: (a) (107.2°, 
-9.388°), (b) (107.4°, -9.205°), (c) (107.6°, -9.648°), (d) (107.7°, -9.411°), (e) (108.3°, 
— 10.02°), (f) (108.2°, —9.75°). The locations of the wave gauges are represented on Figure 
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0.4 
0.3 
0.2 




106 107 108 109 110 

X 

Figure 6. Location of the six numerical wave gauges (indicated by the 
symbol o) superposed with the steady state coseismic bottom displacement 
(only the vertical component in meters is represented here). 



6 along with the static sea bed displacement. Wave gauges are intentionally put in places 
where the largest waves are expected. Synthetic wave gauge records are presented in Figure 
7. We consider the following four scenarios: 

• trigonometric, complete (blue solid line) 

• trigonometric, incomplete (black dashed line) 

• exponential, complete (blue dash-dotted line) 

• exponential, incomplete (black dotted line) 

The importance of the nonlinear effects have already been addressed in several previous 
studies, [52, 76]. In the framework of the Weakly Nonlinear approach we studied this 
question in a recent companion paper [24] using only the vertical deformation. A fairly good 
agreement has been observed with the Cauchy-Poisson (linear, fully dispersive) formulation 
in accordance with preceding results, cf. [52, 76]. Consequently, in the present study we 
decided to focus on the further comparison between complete/incomplete approaches and 
exponential/trigonometric scenarios, [35, 22, 64]. These results are discussed hereafter. 

We note that the trigonometric scenario leads in general to slightly larger amplitudes 
than the exponential bottom motion. It is not surprising since the exponential scenario 
prescribes smoother and less rapid change in the bottom for the same rise time parameter 
value (see Figure 2). Later we will consider the trigonometric scenario unless otherwise 
noted. 

Then, it can be seen that in most cases horizontal displacements lead to an increase in 
the wave amplitude but not always as it can be observed in Figure 7(f) (the last wave 
gauge located at (108.2°, —9.75°)). We investigated more thoroughly this question. Fig- 
ure 8 shows the relative difference between free surface elevations computed according to 
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Gauge 1 




Gauge 2 



Complete, trig. 

Incomplete, trig. 

■ - 1 - 1 Complete, exp. 
Incomplete, exp. 




200 250 



(a) Gauge at (107.2°, -9.388°) 

Gauge 3 



Complete, trig. 
Incomplete, trig. 
Complete, exp. 
Incomplete, exp. 




(b) Gauge at (107.4°, -9.205°) 

Gauge 4 



Complete, trig. 

Incomplete, trig. 

■ - 1 - 1 Complete, exp. 
Incomplete, exp. 




(c) Gauge at (107.6°, -9.648°) 

Gauge 5 



Complete, trig. 
Incomplete, trig. 
Complete, exp. 
Incomplete, exp. 




(d) Gauge at (107.7°, -9.411°) 

Gauge 6 




(e) Gauge at (108.3°, -10.02°) 



(f) Gauge at (108.2°, -9.75°) 



Figure 7. Free surface elevations computed numerically at six wave gauges 
located approximately in local extrema of the static bottom displacement. 
The vertical axis is represented in meters and time is given in seconds. The 
black lines correspond to the wave generated only by the vertical displace- 
ments (incomplete generation) while blue lines take also into account the 
horizontal displacements contribution (complete scenario). 
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-10 




Figure 8. Relative difference between the free surface elevation at t e = 220 
seconds computed according to the complete r) c (x, t e ) (with horizontal dis- 
placements) and incomplete f]i(x,t e ) (only vertical component) tsunami gen- 
eration scenarios. The bottom moves according to the trigonometric scenario 
(2.1). The vertical scale is given in percents of the maximum amplitude of 
the incomplete scenario — max^(f, t e ). 

X 



complete rj c (x,t e ) and incomplete f]i(x,t e ) scenarios. More precisely, we plot the following 
quantity (for the trigonometric scenario): 

= Vc (x,t e ) - v t (x t e ) ^ ^ = g 
max \r)i(x,t e )\ 

X 

Time t e = 220 s has been chosen because at that moment the bottom has been stabilized 
and the waves enter into the free propagation regime. In Figure 8 the grey color corresponds 
to the zero value of the difference d(x), while the white color shows regions where the wave 
is amplified by horizontal displacements by approximately 10%. On the contrary, black 
zones show an attenuation effect of horizontal sea bed motion (about —5%). Recall that 
all values are given in terms of the maximum amplitude max \r)i(x, t e )\ percentage of the 

X 

incomplete generation approach. Some connection with the results presented in Figure 5 
can be noticed. 

Finally, we study the evolution of kinetic, potential and total energies during the tsunami 
generation process described in Section 2.2.1. Specifically we are interested in quantifying 
the contribution of horizontal displacements into tsunami energy balance. Figure 9 shows 
the evolution of potential (9(a)) and kinetic (9(b)) energies for four cases already mentioned 
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Potential energy xlO 12 Kinetic energy 




(a) Potential energy 11(f), J (b) Kinetic energy K(t), J 



Figure 9. Energy evolution during our simulations in the complete (blue 
solid line) and incomplete (black dotted line) scenarios. Note that scales are 
different on the left and right images. The time t is given in seconds. 

above. Here again blue lines refer to complete generation scenarios while black lines - 
to vertical displacements only. Consecutive peaks in the kinetic energy come from the 
activation of new fault segments in accordance with the rupture propagation along the fault. 
The energy curves vary in an analogous way with the tide gauges records (see Figure 7). 
Namely, the trigonometric scenario leads to slightly higher energies than the exponential 
one. However, once the bottom motion stops, this difference becomes negligible. The 
inclusion of horizontal displacements has a much more visible effect with higher energies. 

We performed also the same simulation but using the finite fault solution obtained by 
the Caltech team, [70], and the trigonometric scenario. The results concerning the kinetic 
and potential energies are presented in Figure 10. It is interesting to note that the potential 
energy evolution in the Caltech scenario appears to be smoother especially after t = 150 
s. However, the peaks corresponding to subfaults activation times are equally present in 
kinetic energies. The magnitudes of kinetic and potential energies predicted according 
to the Caltech inversion are approximatively 10 times higher than corresponding USGS 
results. 

We have to underline that inversions performed by the two different finite fault algo- 
rithms lead to very different results. Our method is operational with both of them. It 
is particularly interesting to compare the total energies evolution predicted by USGS and 
Caltech inversions. These results are presented in Figure 11. The Caltech version of the 
rupturing process generates a tsunami wave with much higher energy (computed after the 
end of the bottom motion when a tsunami enters in the propagation regime): 

• USGS, complete trigonometric: E c = 2.16 x 10 12 J, 

• USGS, complete exponential: E c = 2.11 x 10 12 J, 

• USGS, incomplete trigonometric: E{ = 1.95 x 10 12 J, 
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(a) Potential energy 11(f), J (b) Kinetic energy K(t), J 

Figure 10. Energy evolution during our simulations in the complete (blue 
solid line) and incomplete (black dotted line) scenarios using the finite fault 
inversion by Caltech Tectonics Observatory [70]. Note that the scales are 
different on the left and right images. The time t is given in seconds. 



• USGS, incomplete exponential: Ei = 1.90 x 10 12 J, 

• Caltech, complete trigonometric: E c = 2.28 x 10 13 J, 

• Caltech, incomplete trigonometric: E c = 1.83 x 10 13 J. 

This result is not surprising since according to the USGS solution the Java 2006 event 
magnitude is M w = 7.7 and according to the Caltech team, M w = 7.9, which means a huge 
difference since the scale is logarithmic. 

Remark 7. We would like to note an interesting property. For linear waves it can be rigor- 
ously shown the exact equipartition property between the kinetic and potential energies [96]. 
When the nonlinearities are included, this equidistribution property is only approximate. 
One can observe that at the final time in our simulations the kinetic and potential energies 
are already of the same order of magnitude regardless the employed bottom motion (see 
Figures 9 & 10). Moreover, both curves K(t) and H(t) continue to tend to the equilibrium 
state according to the theoretical predictions of the water wave theory, [96] . 

Despite some local attenuation effects of horizontal displacements on the free surface 
elevation, the complete generation scenario produces a tsunami wave with more important 
energy content. More precisely, our computations show that the horizontal displacements 
contribute about 10% into the total tsunami energy balance in the USGS scenario (this 
value is consistent with our previous results concerning the differences in wave amplitudes 
in Figure 8). This result is even more flagrant for the Caltech version which ascribes 24% 
of the energy to horizontal displacements. The free surface amplitudes in this case should 
differ as well by the same order of magnitude. As we already noted, the difference between 
trigonometric and exponential scenarios is negligible. 
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Total energy 



Total energy 





(a) USGS version 



(b) Caltech version 



Figure 1 1 . The total energy evolution predicted according to two different 
versions of the finite fault inversion. Note the different vertical scales on left 
(4.5 x 10 12 ) and right (2.5 x 10 13 ) images. The time t is given in seconds. 

Up to now the contribution of horizontal displacements has been quantified in terms 
of the bottom deformation (Figure 5), free surface elevation at a particular moment of 
time (t = 220 s, cf. Figure 8) and finally in terms of the wave energy (Figures 9 - 
11). However, the wave energy cannot be directly measured in practice. Consequently, 
we continue to quantify the differences between complete and incomplete approaches in 
terms of the waveforms. Namely, we compute also the following relative measures of the 
waveforms difference which act on the whole computed free surface: 

'IVc&t) -r}i(x,t)\\L 2 (M?) , ^ \\Vc(x,t) -r]i(x,t)\\ Lo 



d 2 (t) :-- 



d<x>(t) :-- 



(3.1; 



||fZi(£>*)||za(R2) \\Vi&t)\\ Lo 

The simulation results are presented on Figure 12. One can see that both measures grow up 
to 15% and then oscillate around this level, at least during the simulation time. This result 
is in agreement with previous measurements of the horizontal displacements contribution 
based on the wave energy and the bottom deformation which were of the order of 10%. We 
note also that the curve ^(t) has a more regular behaviour than doo(t) since it is based on 
integral characteristics of the difference while the latter focusses on the characteristics of 
the extreme values. 

It is also interesting to compare the tsunami energy with the energy of the underlying 
seismic event. The USGS Energy and Broadband solution indicates that the radiated 
seismic energy of the Java 2006 earthquake is equal to 3.2 x 10 14 J. Hence, according to 
our computations with the USGS complete generation approach and the trigonometric 
scenario, about 0.68% of the seismic energy was transmitted to the tsunami wave (this 
portion raises to 7.14% for the Caltech version). The ratio of the total tsunami and 
seismic radiated energies can be used as a measure of the tsunami generation efficiency 
of a specific earthquake. The seismic radiated energy may not be the most appropriate 
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Figure 12. Relative difference between free surface elevations computed 
according to the complete and incomplete scenarios. For more details see 
equation (3.1). The horizontal axis is given in seconds, while the vertical 
scale is a percentage. The blue solid line ( — ) corresponds to the measure 
d 2 (t). The black dashed line ( ) represents ^(t). 

parameter to use, but it has an advantage to be relatively robust and easily observable in 
contrast to the more relevant seismic fracture energy, [90, 51]. 

This parameter can also be estimated for the Great Sumatra-Andaman earthquake of 
December 26, 2004. The total energy of the Great Indian Ocean tsunami 2004 was esti- 
mated by T. Lay et al. (2005), [57], to be equal to 4.2 x 10 15 J. According to the same 
USGS Energy and Broadband solution the radiated seismic energy was equal to 1.1 x 10 17 
J. Thus, the tsunamigenic efficiency of the Great Sumatra-Andaman earthquake is 3.81% 
which is bigger than that of Java 2006 event but remains in the same order of magnitude 
around 1%. It is possible that we do not still take into account all important factors in 
Java 2006 tsunami genesis. More examples of transmission of seismic energy to tsunami 
energy can be found e.g. in [58]. 

4. Discussion and conclusions 

In the present study the process of tsunami generation is further investigated. The 
current tsunamigenesis model relies on a combination of the finite fault solution, [45], and a 
recently proposed Weakly Nonlinear (WN) solver for the water wave problem with moving 
bottom, [24]. Consequently, in our model we incorporate recent advances in seismology 
and computational hydrodynamics. 

This study is focused on the role of the horizontal displacements in the real world tsunami 
genesis process. By our intuition we know that a horizontal motion of the flat bottom will 
not cause any significant disturbance on the free surface. However, the real bathymetry 
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is far from being flat. The question which arises naturally is how to quantify the effect of 
horizontal co-seismic displacements during real world events. 

The primary goal of this study was to propose relatively simple, efficient and accurate 
procedures to model tsunami generation process in realistic environments. Special emphasis 
was payed to the role of horizontal displacements which should be also taken into account. 
Thus, the kinematics of horizontal co-seismic displacements were reconstructed and their 
effect on free surface motion was quantified. The evolution of kinetic and potential energies 
were also investigated in our study. In the case of July 17, 2006 Java event our simulations 
indicate that 10% of the energy input can be seemingly ascribed to effects of the horizontal 
bottom motion. This portion increases considerably if we switch to the scenario proposed 
by Caltech Tectonics Observatory [70]. 

The results presented in this study do not still explain the reasons for the extreme run-up 
values caused by July 17, 2006 Java tsunami, [26]. At least we hope that the proposed 
methodology illustrated on this important real world event will be proved helpful in fu- 
ture studies. In our opinion a successful theory should incorporate also other generation 
mechanisms such as local landslides/slumps which are subject to large uncertainties at the 
current stage of our understanding. Until now there are no detailed images of the seabed 
which could support or disprove this assumption. In this study we succeeded to quantify 
the significance of horizontal displacements for the tsunami generation. 

Appendix A. Applications to some recent tsunami events 

In this Appendix we apply the techniques described in this manuscript to two recent 
significant events. We compute the energy transmission from the corresponding seismic 
event to the resulting tsunami wave with and without horizontal displacements. Through- 
out this Appendix we use the trigonometric scenario and the finite fault solutions produced 
by USGS. 

Remark 8. In later versions of the finite fault solution the activation and rise time, as 
well as the seismic moment are specified for each sub-fault. Consequently, below we use 
this information to produce more accurate bottom kinematics. It allows us to compute also 
the seismic moment rate function 5 . 

A.l. Mentawai 2010 tsunami. In October 25, 2010 a small portion of the subduction 
zone seaward of the Mentawai islands was ruptured by an earthquake of magnitude M w = 
7.7, [56]. This earthquake generated a tsunami wave with run-up values ranging from 3 to 
9 meters (with even larger values at some places). On Pagai islands this tsunami caused 
more than 400 victims. This earthquake is characterized by 10° dip angle, a slow rupture 
velocity (~ 1.5 km/s) which propagated during about 100 s over 100 km long source region. 
For our simulation we used the finite fault solution produced by USGS [38]. On Figure 
13 we show the evolution of the total tsunami energy to be compared with the seismic 
moment release rate function. 

5 The total released seismic moment is given by the integral of this function over the rupture duration 
(time). 
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(a) Total energy 



(b) Moment rate 



Figure 13. Total tsunami energy and seismic moment rate function for the 
Mentawai 2010 event. The time t is given in seconds. 



Our computations give the following estimations: 

• Complete scenario: E c = 1.11 x 10 12 J, 

• Incomplete scenario: Ei = 1.06 x 10 12 J. 

Consequently, for this event only 4.5% of the total energy is due to horizontal displacements. 
According to the USGS energy and broadband solution the radiated seismic energy is E s = 
1.4 x 10 15 J. The tsunami energy constitutes 0.08% of the radiated seismic energy. This 
surprising result can be explained by the fact that a big portion of co-seismic displacements 
occured on the land which did not fortunately contribute to the tsunami generation process. 

A. 2. March 11, 2011 Japan earthquake and tsunami. This tsunami event was caused 
by an undersea megathrust earthquake of magnitude M w = 9.0 off the coast of Japan. Offi- 
cially this earthquake was named the Great East Japan Earthquake. The Japanese National 
Police Agency has confirmed more than 13,000 deaths caused both by the earthquake and 
especially by the tsunami. Entire towns were devastated. The local infrastructure includ- 
ing the Fukushima Nuclear Power Plant were heavily affected with consequences which are 
widely known. 

We use again the finite fault inversion performed at the USGS [39]. The tsunami total 
energy evolution along with the moment rate function are represented on Figure 14. At 
the end of the simulation we obtain the following result: 

• Complete scenario: E c = 1.64 x 10 15 J, 

• Incomplete scenario: Ei = 1.50 x 10 15 J. 

Henceforth, the contribution of horizontal displacements is estimated to be 9.2%. The 
radiated seismic energy is E s = 1.9 x 10 17 J according to the USGS energy and broadband 
solution. The total tsunami wave energy represents only 0.87% of the radiated seismic 



CONTRIBUTION OF HORIZONTAL DISPLACEMENTS INTO TSUNAMI GENERATION 



25 




t t, s 

(a) Total energy (b) Moment rate 

Figure 14. Total tsunami energy and seismic moment rate function for the 
Great East Japan Earthquake. The time t is given in seconds. 

energy. This value is much lower than the corresponding value for the Great Sumatra- 
Andaman earthquake of December 26, 2004. 
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